Having established the relationship between NBA performance and betting markets, we now shift our focus to the financial behavior of sports betting stocks themselves. In this chapter, these firms are treated as an independent financial ecosystem shaped by sports events, regulatory changes, and broader market forces. Financial returns exhibit volatility clustering, where periods of high or low volatility tend to persist. Because ARIMA models assume constant variance, they are not well suited for data with this property. GARCH models address this limitation by explicitly capturing time-varying volatility. The key insight is that while returns may be unpredictable due to market efficiency, volatility often follows predictable patterns that are essential for risk management, options pricing, and portfolio optimization.
Price Series: Non-stationary with clear trends. DKNG went public via SPAC merger in 2020, experienced a pandemic boom, then declined as competition intensified.
Returns Series: Appears stationary around zero mean, but exhibits clear volatility clustering. Periods of calm (2021-early 2022) alternate with turbulence (late 2022, 2024).
Code
# Calculate rolling volatility (20-day window)dkng <- dkng %>%arrange(Date) %>%mutate(Roll_Vol_20 = zoo::rollapply(Return_Pct, width =20, FUN = sd, fill =NA, align ="right"),Abs_Return =abs(Return_Pct) )p3 <-ggplot(dkng, aes(x = Date, y = Abs_Return)) +geom_line(color ="purple", alpha =0.6) +labs(title ="Absolute Returns (Proxy for Volatility)",x ="Date", y ="|Return| (%)" ) +theme_minimal()p4 <-ggplot(dkng, aes(x = Date, y = Roll_Vol_20)) +geom_line(color ="darkgreen", linewidth =0.8) +labs(title ="20-Day Rolling Volatility",x ="Date", y ="Volatility (% std dev)" ) +theme_minimal()p3 / p4
Volatility Clustering Evidence:
The rolling volatility plot clearly shows time-varying variance. Volatility is not constant it rises sharply during market stress and falls during stable periods. This visual evidence motivates GARCH modeling.
Code
# ACF/PACF of returnspar(mfrow =c(1, 2))acf(dkng$Return_Pct, lag.max =30, main ="ACF: DKNG Returns")pacf(dkng$Return_Pct, lag.max =30, main ="PACF: DKNG Returns")
Code
par(mfrow =c(1, 1))
Code
# ACF/PACF of squared returns (test for ARCH effects)par(mfrow =c(1, 2))acf(dkng$Return_Pct^2, lag.max =30, main ="ACF: Squared Returns (DKNG)")pacf(dkng$Return_Pct^2, lag.max =30, main ="PACF: Squared Returns (DKNG)")
Code
par(mfrow =c(1, 1))
Key Findings:
Returns ACF/PACF: Minimal autocorrelation at most lags (consistent with weak-form market efficiency). Returns themselves are largely unpredictable.
Squared Returns ACF/PACF: Strong autocorrelation at multiple lags. This is the signature of ARCH effects while returns are unpredictable, volatility is highly predictable and persistent.
Code
# ADF test on returnsadf_dkng <-adf.test(dkng$Return_Pct)cat("=== Augmented Dickey-Fuller Test ===\n")
=== Augmented Dickey-Fuller Test ===
Code
cat("Null Hypothesis: Series has a unit root (non-stationary)\n\n")
Null Hypothesis: Series has a unit root (non-stationary)
cat("Conclusion:", ifelse(adf_dkng$p.value <0.05,"Reject H0 - Returns are STATIONARY","Fail to reject H0 - Returns are NON-STATIONARY"), "\n")
Conclusion: Reject H0 - Returns are STATIONARY
Stationarity Conclusion: The ADF test confirms that DKNG returns are stationary (p < 0.01). This validates our use of returns rather than prices and confirms that we don’t need differencing in the mean equation.
p5 <-ggplot(penn, aes(x = Date, y = Adj_Close)) +geom_line(color ="#006bb6", linewidth =0.8) +labs(title ="Penn Entertainment (PENN) Stock Price",subtitle ="Daily Adjusted Close Price",x ="Date", y ="Price ($)" ) +theme_minimal()p6 <-ggplot(penn, aes(x = Date, y = Return_Pct)) +geom_line(color ="darkred", linewidth =0.5, alpha =0.7) +geom_hline(yintercept =0, linetype ="dashed", color ="black") +labs(title ="Penn Entertainment (PENN) Daily Returns",subtitle ="Log returns showing volatility clustering",x ="Date", y ="Return (%)" ) +theme_minimal()p5 / p6
Observations:
Price Series: Even more volatile than DKNG. PENN’s business model shifted from Barstool Sports (2020-2023) to ESPN BET (2023-present), creating structural breaks and extreme volatility.
Returns Series: Strong volatility clustering with several extreme spikes. This suggests PENN may exhibit even stronger GARCH effects than DKNG.
Code
penn <- penn %>%arrange(Date) %>%mutate(Roll_Vol_20 = zoo::rollapply(Return_Pct, width =20, FUN = sd, fill =NA, align ="right"),Abs_Return =abs(Return_Pct) )p7 <-ggplot(penn, aes(x = Date, y = Abs_Return)) +geom_line(color ="purple", alpha =0.6) +labs(title ="Absolute Returns (Proxy for Volatility)",x ="Date", y ="|Return| (%)" ) +theme_minimal()p8 <-ggplot(penn, aes(x = Date, y = Roll_Vol_20)) +geom_line(color ="darkgreen", linewidth =0.8) +labs(title ="20-Day Rolling Volatility",x ="Date", y ="Volatility (% std dev)" ) +theme_minimal()p7 / p8
Code
par(mfrow =c(1, 2))acf(penn$Return_Pct, lag.max =30, main ="ACF: PENN Returns")pacf(penn$Return_Pct, lag.max =30, main ="PACF: PENN Returns")
Code
par(mfrow =c(1, 1))
Code
par(mfrow =c(1, 2))acf(penn$Return_Pct^2, lag.max =30, main ="ACF: Squared Returns (PENN)")pacf(penn$Return_Pct^2, lag.max =30, main ="PACF: Squared Returns (PENN)")
Code
par(mfrow =c(1, 1))
Key Findings:
Similar to DKNG: returns show weak autocorrelation, but squared returns exhibit strong autocorrelation, indicating ARCH effects. The PACF of squared returns shows significant spikes at early lags, suggesting GARCH(1,1) may be appropriate.
Code
adf_penn <-adf.test(penn$Return_Pct)cat("=== Augmented Dickey-Fuller Test ===\n")
# Use best mean model from AICbest_mean_idx <-which.min(mean_aic$AIC)best_mean_model <- mean_models[[best_mean_idx]]best_mean_name <-names(mean_models)[best_mean_idx]cat("\nSelected Mean Model:", best_mean_name, "\n\n")
Selected Mean Model: ARMA00
Code
# Extract residualsresid_mean <-residuals(best_mean_model)# Plot standardized residuals and squared residualspar(mfrow =c(2, 2))plot(resid_mean, type ="l", main ="Standardized Residuals", ylab ="Residuals")abline(h =0, col ="red", lty =2)plot(resid_mean^2, type ="l", main ="Squared Residuals", ylab ="Residuals")acf(resid_mean^2, lag.max =30, main ="ACF: Squared Residuals")pacf(resid_mean^2, lag.max =30, main ="PACF: Squared Residuals")
penn_returns <- penn$Return_Pct# Fit candidate mean modelsmean_models_penn <-list(ARMA00 =arima(penn_returns, order =c(0, 0, 0)),ARMA10 =arima(penn_returns, order =c(1, 0, 0)),ARMA01 =arima(penn_returns, order =c(0, 0, 1)),ARMA11 =arima(penn_returns, order =c(1, 0, 1)))mean_aic_penn <-data.frame(Model =names(mean_models_penn),AIC =sapply(mean_models_penn, AIC),BIC =sapply(mean_models_penn, BIC))cat("=== PENN Mean Equation Comparison ===\n\n")
=== PENN Mean Equation Comparison ===
Code
kable(mean_aic_penn,format ="html",digits =2,caption ="Model Comparison: Mean Equation Only (PENN)") %>%kable_styling(full_width =FALSE,bootstrap_options =c("striped", "hover", "condensed", "responsive") ) %>%row_spec(which.min(mean_aic_penn$AIC), bold =TRUE, color ="white", background ="#006bb6")
garch_models_penn <-list()garch_results_penn <-data.frame()for (spec in model_specs) {# Fit model using fGarch fit <-tryCatch( {garchFit(formula = spec$formula, data = penn_returns, trace =FALSE) },error =function(e) {cat("Warning: Failed to fit", spec$name, "for PENN\n")NULL } )if (!is.null(fit)) { garch_models_penn[[spec$name]] <- fit# Extract information criteria aic_val <- fit@fit$ics["AIC"] bic_val <- fit@fit$ics["BIC"] loglik_val <- fit@fit$llh garch_results_penn <-rbind(garch_results_penn, data.frame(Model = spec$name,AIC = aic_val,BIC = bic_val,Log_Likelihood = loglik_val )) }}cat("=== PENN GARCH Model Comparison ===\n\n")
=== PENN GARCH Model Comparison ===
Model Comparison: ARMA+GARCH(1,1) Models for PENN
Model
AIC
BIC
Log_Likelihood
ARMA(0,0)+GARCH(1,1)
5.6132
5.6295
3523.881
ARMA(1,0)+GARCH(1,1)
5.6141
5.6345
3523.455
ARMA(0,1)+GARCH(1,1)
5.6141
5.6345
3523.469
ARMA(1,1)+GARCH(1,1)
5.6156
5.6402
3523.433
Code
best_garch_penn <- garch_models_penn[[which.min(garch_results_penn$AIC)]]best_name_penn <- garch_results_penn$Model[which.min(garch_results_penn$AIC)]cat("\n*** BEST MODEL:", best_name_penn, "***\n\n")
*** BEST MODEL: ARMA(0,0)+GARCH(1,1) ***
Code
cat("Model Coefficients:\n")
Model Coefficients:
Code
print(coef(best_garch_penn))
mu omega alpha1 beta1
-0.03411880 0.86283109 0.09135361 0.86147898
Code
std_resid_penn <-residuals(best_garch_penn, standardize =TRUE)par(mfrow =c(3, 2))plot(std_resid_penn, type ="l", main ="Standardized Residuals (PENN)")abline(h =0, col ="red", lty =2)plot(std_resid_penn^2, type ="l", main ="Squared Std. Residuals (PENN)")acf(std_resid_penn, lag.max =20, main ="ACF: Std. Residuals")acf(std_resid_penn^2, lag.max =20, main ="ACF: Squared Std. Residuals")qqnorm(std_resid_penn, main ="Q-Q Plot")qqline(std_resid_penn, col ="red")hist(std_resid_penn, breaks =50, probability =TRUE, main ="Distribution")curve(dnorm(x), add =TRUE, col ="red", lwd =2)
Code
par(mfrow =c(1, 1))# Ljung-Box testslb_penn <-Box.test(std_resid_penn, lag =10, type ="Ljung-Box")lb_penn_sq <-Box.test(std_resid_penn^2, lag =10, type ="Ljung-Box")cat("\nLjung-Box (Std. Residuals): p =", round(lb_penn$p.value, 4), "\n")
Ljung-Box (Std. Residuals): p = 0.438
Code
cat("Ljung-Box (Squared): p =", round(lb_penn_sq$p.value, 4), "\n")
Interpretation: The ACF/PACF of squared residuals show significant autocorrelation at multiple lags, and the ARCH LM test confirms the presence of ARCH effects. This validates the need for GARCH modeling to capture CZR’s time-varying volatility.
Code
garch_models_czr <-list()garch_results_czr <-data.frame()for (spec in model_specs) {# Fit model using fGarch fit <-tryCatch( {garchFit(formula = spec$formula, data = czr_returns, trace =FALSE) },error =function(e) {cat("Warning: Failed to fit", spec$name, "for CZR\n")NULL } )if (!is.null(fit)) { garch_models_czr[[spec$name]] <- fit# Extract information criteria aic_val <- fit@fit$ics["AIC"] bic_val <- fit@fit$ics["BIC"] loglik_val <- fit@fit$llh garch_results_czr <-rbind(garch_results_czr, data.frame(Model = spec$name,AIC = aic_val,BIC = bic_val,Log_Likelihood = loglik_val )) }}cat("=== CZR GARCH Model Comparison ===\n\n")
=== CZR GARCH Model Comparison ===
Model Comparison: ARMA+GARCH(1,1) Models for CZR
Model
AIC
BIC
Log_Likelihood
ARMA(0,0)+GARCH(1,1)
5.3783
5.3926
3967.845
ARMA(1,0)+GARCH(1,1)
5.3784
5.3963
3966.944
ARMA(0,1)+GARCH(1,1)
5.3784
5.3963
3966.913
ARMA(1,1)+GARCH(1,1)
5.3796
5.4011
3966.849
Code
best_garch_czr <- garch_models_czr[[which.min(garch_results_czr$AIC)]]best_name_czr <- garch_results_czr$Model[which.min(garch_results_czr$AIC)]cat("\n*** BEST MODEL:", best_name_czr, "***\n\n")
*** BEST MODEL: ARMA(0,0)+GARCH(1,1) ***
Code
cat("Model Coefficients:\n")
Model Coefficients:
Code
print(coef(best_garch_czr))
mu omega alpha1 beta1
-0.01029034 0.66993148 0.11243921 0.84204206
Model Selection Rationale: The AIC-selected model balances fit and parsimony. CZR exhibits volatility patterns driven by both established gaming operations and newer sports betting ventures.
Code
std_resid_czr <-residuals(best_garch_czr, standardize =TRUE)par(mfrow =c(3, 2))plot(std_resid_czr, type ="l", main ="Standardized Residuals (CZR)", ylab ="Std. Residuals")abline(h =0, col ="red", lty =2)abline(h =c(-2, 2), col ="blue", lty =2)plot(std_resid_czr^2, type ="l", main ="Squared Std. Residuals (CZR)", ylab ="Residuals²")acf(std_resid_czr, lag.max =20, main ="ACF: Std. Residuals")acf(std_resid_czr^2, lag.max =20, main ="ACF: Squared Std. Residuals")qqnorm(std_resid_czr, main ="Q-Q Plot (CZR)")qqline(std_resid_czr, col ="red")hist(std_resid_czr, breaks =50, probability =TRUE, main ="Distribution (CZR)", xlab ="Std. Residuals")curve(dnorm(x), add =TRUE, col ="red", lwd =2)
Code
par(mfrow =c(1, 1))# Ljung-Box testslb_czr <-Box.test(std_resid_czr, lag =10, type ="Ljung-Box")lb_czr_sq <-Box.test(std_resid_czr^2, lag =10, type ="Ljung-Box")cat("\n=== Ljung-Box Tests (CZR) ===\n\n")
=== Ljung-Box Tests (CZR) ===
Code
cat("Standardized Residuals: p =", round(lb_czr$p.value, 4), "\n")
Volatility Interpretation: CZR’s conditional volatility shows periods of elevated risk during major corporate events. The GARCH model successfully tracks these time-varying patterns, providing valuable risk forecasts for portfolio managers and options traders.
Interpretation: Despite MGM’s diversified business model, the stock still exhibits ARCH effects, though potentially milder than pure-play sports betting stocks. The ARCH LM test confirms time-varying volatility requiring GARCH specification.
Code
garch_models_mgm <-list()garch_results_mgm <-data.frame()for (spec in model_specs) {# Fit model using fGarch fit <-tryCatch( {garchFit(formula = spec$formula, data = mgm_returns, trace =FALSE) },error =function(e) {cat("Warning: Failed to fit", spec$name, "for MGM\n")NULL } )if (!is.null(fit)) { garch_models_mgm[[spec$name]] <- fit# Extract information criteria aic_val <- fit@fit$ics["AIC"] bic_val <- fit@fit$ics["BIC"] loglik_val <- fit@fit$llh garch_results_mgm <-rbind(garch_results_mgm, data.frame(Model = spec$name,AIC = aic_val,BIC = bic_val,Log_Likelihood = loglik_val )) }}cat("=== MGM GARCH Model Comparison ===\n\n")
=== MGM GARCH Model Comparison ===
Model Comparison: ARMA+GARCH(1,1) Models for MGM
Model
AIC
BIC
Log_Likelihood
ARMA(0,0)+GARCH(1,1)
4.9152
4.9296
3625.892
ARMA(1,0)+GARCH(1,1)
4.9165
4.9344
3625.799
ARMA(0,1)+GARCH(1,1)
4.9165
4.9344
3625.800
ARMA(1,1)+GARCH(1,1)
4.9142
4.9357
3623.114
Code
best_garch_mgm <- garch_models_mgm[[which.min(garch_results_mgm$AIC)]]best_name_mgm <- garch_results_mgm$Model[which.min(garch_results_mgm$AIC)]cat("\n*** BEST MODEL:", best_name_mgm, "***\n\n")
Model Selection Rationale: The selected GARCH(1,1) specification for MGM captures volatility persistence while maintaining parsimony. MGM’s diversified operations (Las Vegas Strip properties, regional casinos, BetMGM, entertainment venues) create a more stable volatility profile compared to pure-play operators, which should be reflected in the model coefficients.
Code
std_resid_mgm <-residuals(best_garch_mgm, standardize =TRUE)par(mfrow =c(3, 2))plot(std_resid_mgm, type ="l", main ="Standardized Residuals (MGM)", ylab ="Std. Residuals")abline(h =0, col ="red", lty =2)abline(h =c(-2, 2), col ="blue", lty =2)plot(std_resid_mgm^2, type ="l", main ="Squared Std. Residuals (MGM)", ylab ="Residuals²")acf(std_resid_mgm, lag.max =20, main ="ACF: Std. Residuals")acf(std_resid_mgm^2, lag.max =20, main ="ACF: Squared Std. Residuals")qqnorm(std_resid_mgm, main ="Q-Q Plot (MGM)")qqline(std_resid_mgm, col ="red")hist(std_resid_mgm, breaks =50, probability =TRUE, main ="Distribution (MGM)", xlab ="Std. Residuals")curve(dnorm(x), add =TRUE, col ="red", lwd =2)
Code
par(mfrow =c(1, 1))# Ljung-Box testslb_mgm <-Box.test(std_resid_mgm, lag =10, type ="Ljung-Box")lb_mgm_sq <-Box.test(std_resid_mgm^2, lag =10, type ="Ljung-Box")cat("\n=== Ljung-Box Tests (MGM) ===\n\n")
=== Ljung-Box Tests (MGM) ===
Code
cat("Standardized Residuals: p =", round(lb_mgm$p.value, 4), "\n")
Diagnostic Summary: MGM’s standardized residuals should show cleaner behavior than PENN or DKNG, reflecting the stabilizing effect of diversified revenue streams. Successful diagnostics (Ljung-Box p > 0.05) indicate the GARCH model adequately captures MGM’s conditional volatility dynamics.
Volatility Interpretation: MGM’s conditional volatility time series shows the lowest average volatility among all four stocks, confirming that diversification provides natural hedging against sector-specific shocks. The GARCH forecasts are particularly valuable for MGM during major events that create temporary volatility spikes.
Unconditional volatility: 2.9379% (lowest among all four stocks)
Business Model Impact: MGM’s lower volatility parameters directly reflect its diversified operations. The smaller \(\alpha_1\) suggests that MGM’s stock is less sensitive to short-term shocks, while moderate \(\beta_1\) indicates reasonable volatility persistence.
Cross-Stock Comparison
Now that we’ve analyzed individual stocks, let’s compare volatility dynamics across all four sports betting companies.
The financial models above demonstrate that sports betting stocks exhibit classic volatility clustering patterns captured by ARCH/GARCH frameworks. DKNG, PENN, CZR, and MGM show time-varying variance where turbulent periods persist before giving way to calm. These models excel at forecasting conditional volatility—predicting the uncertainty around returns rather than the returns themselves—making them essential for risk management and options pricing.
Yet ARCH/GARCH models share a fundamental assumption with ARIMA and VAR: they impose explicit parametric structure. GARCH(1,1) assumes today’s variance depends on yesterday’s variance and yesterday’s squared return with fixed coefficients. ARIMA(p,d,q) requires us to specify lag orders before fitting. These models demand domain expertise and statistical judgment, making them powerful when assumptions hold but brittle when reality deviates.
The deep learning alternative abandons explicit parameterization in favor of learned representations. Recurrent neural networks (RNNs), Long Short-Term Memory (LSTM), and Gated Recurrent Units (GRUs) automatically discover temporal patterns through gradient descent rather than requiring analysts to pre-specify AR and MA terms. This flexibility enables deep learning models to capture complex nonlinearities, long-range dependencies, and regime shifts that defy simple parametric forms.
The next chapter tests whether this flexibility improves forecasts or simply overfits limited data. With only 45 years of NBA data, do LSTM’s 1,600+ parameters extract genuine signal or memorize noise? The comparison is direct: identical train/test splits, identical evaluation metrics (RMSE, MAE), identical target variable (ORtg). Traditional methods have interpretable coefficients and well-understood statistical properties; deep learning offers representational power but requires careful regularization.
What you’ll see next:
RNN, LSTM, and GRU architectures applied to NBA offensive efficiency forecasting, with explicit comparison to ARIMA baselines
Regularization techniques (L2 penalties, dropout, early stopping) that prevent overfitting on small datasets
Input windowing that transforms annual observations into supervised learning sequences (5-year windows → 1-year ahead)
Performance comparison tables showing whether deep learning adds value beyond traditional methods for sports analytics
The transition from financial to deep learning models reflects a methodological tension in modern forecasting: parametric parsimony versus representational flexibility. GARCH models have 2-3 parameters; LSTMs have thousands. For volatile betting stocks with rich intraday data, that complexity may capture genuine structure. For annual NBA statistics with 45 observations, it may simply overfit.
The critical question: Does basketball’s evolution exhibit patterns too complex for ARIMA’s linear assumptions but too simple for LSTM’s deep architecture? Or do we sit in the sweet spot where both approaches offer complementary insights? The empirical comparison tests whether sports analytics has entered the deep learning era or whether classical time series methods remain the gold standard for data-scarce domains.
Where GARCH asked “How does volatility evolve?”, deep learning asks “What temporal patterns can we learn without assuming functional form?” Both frameworks forecast the future, but from opposite epistemological starting points: theory-driven specification versus data-driven discovery.